


options_fmin=struct('Display','iter','MaxFunEvals',100000,'MaxIter',1000,'GradObj','off','LargeScale','off','UseParallel',false,'TolFun',1e-6);

%Set up data
b = [0.060 0.165 0.068 0.037 0.086 0.024]';
se = [0.030 0.084 0.072 0.053 0.031 0.012]';
V= diag(se.^2);

%Fit cmd model
mu_0=mean(b);
[mu_hat Q_hat exitflag output G H]=fminunc(@(x)(b-x)'*inv(V)*(b-x),mu_0,options_fmin);
	V=inv(H);
	SE=sqrt(diag(V));
	[mu_hat SE]
	1-chi2cdf(Q_hat,length(b)-1)

%Set up data
b = [0.054 0.170 0.092 0.057 0.054]';
se = [0.029 0.068 0.056 0.036 0.028]';
V= diag(se.^2);

%Fit cmd model
mu_0=mean(b);
[mu_hat Q_hat exitflag output G H]=fminunc(@(x)(b-x)'*inv(V)*(b-x),mu_0,options_fmin);
	V=inv(H);
	SE=sqrt(diag(V));
	[mu_hat SE]
	1-chi2cdf(Q_hat,length(b)-1)
 

%Set up data
%b = [0.052 0.170 0.120]';
%se = [0.029 0.068 0.031]';
%V= diag(se.^2);

%Fit cmd model
%mu_0=mean(b);
%[mu_hat Q_hat exitflag output G H]=fminunc(@(x)(b-x)'*inv(V)*(b-x),mu_0,options_fmin);
%	V=inv(H);
%	SE=sqrt(diag(V));
%	[mu_hat SE]
%	1-chi2cdf(Q_hat,length(b)-1)